指数関数のパラメータの推定
a2,をシフト
\(\Large \displaystyle y_i = a_0 \ Exp (- a_1 x_i) + a_2 \)
a0 | 10.5226 | 10.4345 | 10.3454 | 10.2552 | 10.1640 | 10.0718 | 9.9786 | 9.8845 | 9.7894 | |||
a1 | 0.4462 | 0.4577 | 0.4696 | 0.4819 | 0.4946 | 0.5077 | 0.5213 | 0.5355 | 0.5501 | |||
a2 | 0.4839 | 0.5839 | 0.6839 | 0.7839 | 0.8839 | 0.9839 | 1.0839 | 1.1839 | 1.2839 | |||
δ | -0.4 | -0.3 | -0.2 | -0.1 | 0 | 0.1 | 0.2 | 0.3 | 0.4 | |||
i | x | y | \( \hat{y} \) | |||||||||
1 | 0 | 10 | 11.0065 | 11.0185 | 11.0293 | 11.0391 | 11.0479 | 11.0557 | 11.0625 | 11.0684 | 11.0733 | |
2 | 2 | 4 | 4.7943 | 4.7612 | 4.7284 | 4.6960 | 4.6639 | 4.6324 | 4.6015 | 4.5712 | 4.5416 | |
3 | 3 | 2 | 3.2427 | 3.2270 | 3.2128 | 3.2001 | 3.1891 | 3.1799 | 3.1724 | 3.1668 | 3.1632 | |
4 | 4 | 1 | 2.2496 | 2.2563 | 2.2651 | 2.2763 | 2.2897 | 2.3056 | 2.3239 | 2.3447 | 2.3680 | |
5 | 6 | 0.5 | 1.2072 | 1.2534 | 1.3021 | 1.3532 | 1.4067 | 1.4627 | 1.5210 | 1.5817 | 1.6447 | |
6 | 9 | 0.1 | 0.6735 | 0.7535 | 0.8350 | 0.9181 | 1.0025 | 1.0883 | 1.1754 | 1.2637 | 1.3532 | |
S (\(y_i - \hat{y} \)の平方和) | 0.4311 | 0.3554 | 0.2996 | 0.2650 | 0.2531 |
0.2655 | 0.3035 | 0.3687 | 0.4626 | |||
dS (Seとの差分) | 0.1780 | 0.1023 | 0.0464 | 0.0118 | 0 | 0.0123 | 0.0503 | 0.1155 | 0.2094 | |||
・残差平方和
推定値からの残差
\(\Large \displaystyle Se = \sum_{i=1}^{n} \left[ y_i -\hat{a_0} \ Exp(- \hat{a_1} \ x_i) - \hat{a_2} \right]^2 \)
a0をシフトさせたときの,推定値からの残差
\(\Large \displaystyle Se = \sum_{i=1}^{n} \left[ y_i -\hat{a_0} \ Exp(- \hat{a_1} \ x_i) -a_2 \right]^2 \)
であり,a1を,δ,だけシフトさせて,固定し,その際のa0, a2の推定値をソルバーで推定しました.
dS,を見ていただけるとわかるように,推定値,Seが一番小さく,今回は,ほぼ左右対称に増加していることがわかります.
グラフ化すると,
のように,二乗+定数できれいに近似できます.
二次曲線の近似においてもきれいに近似でき,
\(\Large \displaystyle y = 0.2531 + 1.2105 \ \delta^2 \)
ここで,分散値は,
・分散
\(\Large \displaystyle Ve = \frac{1}{n-3} \sum_{i=1}^{n} \left[ y_i -\hat{a_0} \ Exp(- \hat{a_1} \ x_i) \right]^2 = \frac{Se}{n-3} = \frac{0.2531}{6-3} = 0.08438 \)
であり(a0,a1,a2の二つのパラメータが3つあるので,自由度は,n-3),
\(\Large \displaystyle 1.2105 \ \delta^2 = 0.08438 \)
となるδがSEとなるので,
\(\Large \displaystyle \delta^2 = \frac{ 0.08438}{1.2105} =0.06971 \)
\(\Large \displaystyle SE_{a_0} = \sqrt{\delta^2} =\color{red}{0.2640} \)
と推定できます.
・Rによる推定
Rでの近似を行ってみると,
プログラムは,
xx <- c(0,2,3,4,6,9)
yy <- c(11,5,3,2,1.5,1.1)
plot(xx,yy)
fm<-nls(yy~a0*exp(-a1*xx)+a2,start=c(a0=10,a1=0.5,a2=1),trace=TRUE)
summary(fm)
で,結果は,
Parameters: | |||||
Estimate | Std. Error | t value | Pr(>|t|) | ||
a0 | 10.16402 | 0.37708 | 26.954 | 0.000112 | *** |
a1 | 0.49456 | 0.04408 | 11.219 | 0.001518 | ** |
a2 | 0.88393 | 0.26931 | 3.282 | 0.046348 | * |
となり,Kyplotにおいても,
推定値 | 標準誤差(SE) | |
A1 | 10.16401 | 0.377081 |
A2 | 0.494562 | 0.044083 |
A3 | 0.883926 | 0.269308 |
と同じ結果となり,今回の推定値と,ほぼ一致,します.
微妙に異なるのが気になりますが....
「統計解析の初歩」,の「1.2 非線形最小2乗法の基本的な考え方」には,
δのずらした値0.5 を変えると結果は異なり、近似標準誤差の精度が変わる
統計パッケージにより、近似標準誤差の値は幾分異なる
とあります.ここで,”0.5”,がどこから出てきたかはわかりません.そもそも横軸(x軸)の範囲に依存しちゃいますし..
そこで,σをとる値(±での平均値)でどう推定値が変わるかを調べてみました.その結果が,
とどのδでもRなどの推定値より下回っていました....謎です...
次ページに今までのまとめを記します.